%matplotlib inline
Calculations of distributions of proplyd distances and sizes and angles. These translate into distributions of $\beta$.
Data sources¶
We have the list of proplyds from Henney & Arthur (1998) that we used in Fig 12 of Tsamis et al (2013). This has distances and sizes.
We also have the Ricci et al (2008) catalog, which has distances but not sizes.
from astropy.table import Table
import numpy as np
import seaborn as sns
from scipy import stats
Now use seaborn to get nicer graphs.
sns.set_palette("deep", desat=0.6)
sns.set_context("talk", rc={"figure.figsize": (12, 8)})
Henney & Arthur (1998) sample¶
datafile = "Proplyd-Datasets/vicente-vs-HA98.dat"
data = Table.read(datafile, format='ascii')
data
from matplotlib import pyplot as plt
_ = plt.hist(data["d"], bins=30, range=(0, 30), cumulative=True)
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of proplyds with projected distance < D\'')
_ = plt.title('Cumulative distribution of proplyds vs distance: HA98')
m1 = data["d"] <= 7.5
m2 = np.abs(data["d"] - 13.0) <= 5.5
m3 = np.abs(data["d"] - 25.0) <= 6.5
h1, edges = np.histogram(data["r14"][m1], bins=10, range=(0, 30), density=True)
h2, edges = np.histogram(data["r14"][m2], bins=10, range=(0, 30), density=True)
h3, edges = np.histogram(data["r14"][m3], bins=10, range=(0, 30), density=True)
plt.bar(edges[:-1], h1, 2.5, color='r', alpha=0.5, label='close: D\' < 7.5"')
plt.bar(edges[:-1], h2, 2.5, h1, color='g', alpha=0.5, label='medium: 7.5" < D\' < 18.5"')
plt.bar(edges[:-1], h3, 2.5, h1 + h2, color='b', alpha=0.5, label='far: D\' > 18.5"')
_ = plt.xlabel('Proplyd i-front radius, 1e14 cm')
_ = plt.ylabel('Fraction of proplyds in each size bin')
_ = plt.title('Histogram of proplyd sizes: HA98')
_ = plt.legend(loc='upper right')
_ = plt.plot(data["d"], data["r14"], 'o')
_ = plt.xlim(0.0, 30.0)
_ = plt.ylim(0.0, 30.0)
plt.xlabel("Distance, arcsec")
plt.ylabel("Radius, 1e14 cm")
We need to take account of binaries
Now we will try out some of the seaborn tools for visualizing data distributions. See this tutorial.
In the section on the Ricci data it occured to me that we need to mirror the data about zero because the radius and separation are positive-definite.
def mirrored(dataset):
"""Expand a positive-definite dataset (1d sequence) to include its negative mirror set"""
assert np.all(dataset >= 0.0)
return np.r_[dataset, -dataset]
def doubled(dataset):
"""Double up a dataset (1d sequence)"""
return np.r_[dataset, dataset]
m30 = data["d"] <= 30.0
with sns.axes_style("white"):
sns.jointplot(mirrored(data["d"][m30]), mirrored(data["r14"][m30]),
kind="hex", size=9, xlim=(0, 30), ylim=(0, 30))
sns.kdeplot(mirrored(data["d"]), mirrored(data["r14"]),
clip=((0, 40), (0, 30)),
bw=3.0, shade=True, cmap="Reds",
)
None
color = "#774499"
g = sns.JointGrid(mirrored(data["d"][m30]), mirrored(data["r14"][m30]), size=9, xlim=(0, 30), ylim=(0, 30))
# g.plot(sns.regplot, sns.distplot, stats.pearsonr);
g.plot_marginals(sns.distplot, color=color, bins=10, rug=True, kde_kws=dict(bw=5.0, label="KDE"))
g.plot_joint(sns.regplot, color=color, ci=98, order=3, label="Cubic regression fit with 98% confidence bounds")
g.ax_joint.legend(loc="upper left")
g.ax_joint.set_xlabel("Projected separation, $D'$")
g.ax_joint.set_ylabel("Proplyd radius $r_0$, $10^{14}$ cm")
g.ax_marg_x.legend(loc="upper right")
g.ax_marg_y.legend(loc="best")
g.fig.tight_layout()
g.fig.savefig("proplyd-joint-separation-size.pdf");
#g.annotate(stats.pearsonr);
Now we will compare the sizes for different intervals of $D'$
datalist = [data['r14'][m1], data['r14'][m2], data['r14'][m3]]
mirrored_datalist = [mirrored(x) for x in datalist]
namelist = ['close: D\' < 7.5"', 'mid: 7.5" < D\' < 18.5"', 'far: D\' > 18.5"']
The simplest option is boxplot, which shows the median and quartiles. By setting whis=np.inf we make the whiskers extend out to the full range of the data.
sns.boxplot(datalist, names=namelist, whis=np.inf, color='pastel')
plt.ylabel("Proplyd size, 1e14 cm")
plt.xlabel("Projected distance band")
None
Alternatively we can use violinplot, which shows the Kernel Density Estimate (KDE) of the values in each category. By setting inner="points" it shows the actual data points instead of the quartiles.
I am a little uncomfortable of this way of presenting the data, since we are forcing the distance to be categorical, when really it is continuous.
sns.violinplot(mirrored_datalist,
names=namelist,
inner="points", bw=0.3, color="PuBuGn")
plt.ylim(0.0, None)
plt.ylabel("Proplyd size, 1e14 cm")
plt.xlabel("Projected distance band")
None
Ricci et al (2008) sample¶
rdata = Table.read("Proplyd-Datasets/ricci-2008.dat", format='ascii')
rdata.remove_row(0) # First row is rubbish
rdata
Bright proplyds have i in the Type column.
m = rdata['Type'] == 'i'
md = rdata['Type'] == 'd'
mj = rdata['Type'] == 'j'
mr = rdata['Type'] == 'rn'
rdata[m]
print("Number of bright proplyds:", len(rdata[m]))
from astropy import coordinates as coord
import astropy.units as u
Try out the new coordinates interface.
c = coord.SkyCoord(ra="05 35 33.20", dec="-05 16 05.38", unit=("hourangle", "deg"))
c
Add a column with an astropy SkyCoord instance for each object
rdata["coord"] = coord.SkyCoord(ra=rdata["RAJ2000"], dec=rdata["DEJ2000"], unit=("hourangle", "deg"))
Find the coordinates of $\theta^1$ Ori C:
c0 = coord.get_icrs_coordinates("tet01 ori c")
c0
c0.separation(c).arcsec
Add another column with the projected separation in arcsec:
rdata['Dprime'] = [c0.separation(c).arcsec for c in rdata["coord"]]
Make a cumulative histogram to compare with the HA98 sample:
rdata['Dprime'][m].max(), rdata['Dprime'].max()
rmax = 800.0
nbins = rmax//4
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=True, color='b', alpha=0.3, label='Ricci 2008')
_ = plt.hist(data["d"], bins=nbins, range=(0, rmax), cumulative=True, color='r', alpha=0.3, label='Henney & Arthur 1998')
_ = plt.hist(rdata["Dprime"][md], bins=nbins, range=(0, rmax), cumulative=True, color='y', alpha=0.7, label='Dark disks')
_ = plt.hist(rdata["Dprime"][mr], bins=nbins, range=(0, rmax), cumulative=True, color='g', alpha=0.7, label='Reflec nebs')
_ = plt.hist(rdata["Dprime"][mj], bins=nbins, range=(0, rmax), cumulative=True, color='m', alpha=0.7, label='Bare jets')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of proplyds with projected distance < D\'')
_ = plt.title('Cumulative distribution of proplyds vs distance: Both samples')
plt.legend(loc='upper left')
plt.gcf().set_size_inches(18, 8)
plt.gcf().savefig("proplyd-distance-distribution.pdf")
Try it again, but non-cumulative. If the cumulative distribution is linear in $D'$ then the non-cumulative one should be constant.
rmax = 800.0
nbins = rmax//20
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=False, color='b', alpha=0.3, label='Ricci 2008')
_ = plt.hist(data["d"], bins=nbins, range=(0, rmax), cumulative=False, color='r', alpha=0.3, label='Henney & Arthur 1998')
_ = plt.hist(rdata["Dprime"][md], bins=nbins, range=(0, rmax), cumulative=False, color='y', alpha=0.7, label='Dark disks')
_ = plt.hist(rdata["Dprime"][mr], bins=nbins, range=(0, rmax), cumulative=False, color='g', alpha=0.7, label='Reflec nebs')
_ = plt.hist(rdata["Dprime"][mj], bins=nbins, range=(0, rmax), cumulative=False, color='m', alpha=0.7, label='Bare jets')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of proplyds with projected distance < D\'')
_ = plt.title('Distribution of proplyds vs distance: Both samples')
plt.legend(loc='upper right')
plt.gcf().set_size_inches(12, 6)
The ultimate tool for plotting a distribution is distplot. This can show the following if we turn on all the options:
- A rugplot of the individual data points
- A histogram of the distribution
- The KDE of the distribution
- A fit to the distribution of any function from
scipy.stats
sns.distplot(rdata['Dprime'][m],
rug=True,
fit=stats.expon,
kde_kws={"clip": (0.0, 800.0), "kernel": "gau", "label": "Gaussian KDE"},
rug_kws={"alpha": 0.3},
fit_kws={"label": "Exponential fit"},
bins=800//25,
label="Histogram (25 arcsec bins)", axlabel=r"Projected distance, $D'$")
plt.legend()
plt.xlim(0.0, 800.0)
plt.title("Distribution of projected distances for all bright proplyds from Ricci et al. (2008) sample")
None
What I don't like about this is that the KDE is affected by the fact that there are no points with $D' < 0$. That is why it drops towards zero. To fix this, we could reflect the dataset about zero and add a negative mirrored dataset.
Note that stats.laplace is actually a 2-sided explonential. And it doesn't fit very well. On the other hand, the KDE now looks much nicer.
sns.distplot(np.r_[rdata['Dprime'][m], -rdata['Dprime'][m]],
rug=True,
fit=stats.logistic,
kde_kws={"kernel": "gau", "label": "KDE"},
rug_kws={"alpha": 0.3},
fit_kws={"label": "Sech-squared fit"},
bins=2*800//25,
label="Histogram (25 arcsec bins)", axlabel=r"Projected distance, $D'$")
plt.legend()
plt.xlim(0.0, 800.0)
plt.title("Distribution of projected distances for all bright proplyds from Ricci et al. (2008) sample")
None
So now, we will try the same trick on the joint distribution of $D'$ and $r_0$. This is in a previous section.
Zooming in on the central part, which is relevant to the bowshocks.
Dmax = 110.0
Dmaxx = 800.0
mzoom = rdata['Dprime'] < Dmax
dataset = mirrored(np.array(rdata['Dprime'][m & mzoom]))
sns.distplot(dataset,
rug=True,
fit=stats.uniform,
kde_kws={"kernel": "gau", "label": "Gaussian KDE"},
rug_kws={"alpha": 0.3},
fit_kws={"label": "Constant fit"},
hist_kws={"range": (-Dmaxx, Dmaxx)},
bins=2*Dmaxx//10,
label="Histogram (10 arcsec bins)", axlabel=r"Projected distance, $D'$")
plt.legend(loc="upper left")
plt.xlim(0.0, 100.0)
plt.title("Distribution of projected distances for all bright proplyds from Ricci et al. (2008) sample");
So this suggests that at scales $< 100''$ the proplyd distribution is uniform (implying density $\sim D^{-2}$.
Whereas at larger scales it does fall off.
What about the stars?¶
sdata = Table.read("Proplyd-Datasets/Robberto2013/table5.dat", format='ascii')
sdata.colnames
We first need to restrict the table to unique sources, since it contains multiple observations of the same source. We can do this by comparing the first two columns:
m_uniq = sdata['Seq'] == sdata['oncacs']
Then winnowing down the table:
sdata = sdata[m_uniq]
Description of the type column from notes to Table 5:
Source type, from visual inspection: (0) not measured; (1) detected in at least one filter and unresolved; (2) double (companion closer than ≈3 pixels, one entry for both sources); (3) wide double (companion further than ≈3 pixel, one entry for each source); (5) silhouette disk; (6) photoionized (proplyd or with other evidence of photoionization); (7) galaxy; (8) Herbig-Haro.source_types = ['not measured', 'single', 'close double', 'wide double',
'???', 'disk', 'ionized', 'galaxy', 'herbig haro']
count_table = Table(names=['Code', 'Type', 'Number'],
dtype=[int, 'S12', int])
for itype, label in enumerate(source_types):
count_table.add_row([itype, label, np.sum(sdata['type'] == itype)])
count_table
The ionized type is "proplyd or with other evidence of photoionization" and there are 184 of them, compared with 178 in the Ricci catalog. This is close enough for me.
Next job is to winnow out types 0, 7, and 8, which we don't want.
sdata = sdata[(sdata['type'] > 0) & (sdata['type'] < 7)]
And now create a mask to select out the ionized sources (mainly proplyds).
m_rob_ionized = sdata['type'] == 6
Now we combine the 6 columns used for coordinates into a single RA-Dec string:
sdata["coordstring"] = ["{} {} {} -{} {} {}".format(*args)
for args in zip(*[sdata[col] for col in ["RAh", "RAm", "RAs", "DEd", "DEm", "DEs"]])]
sdata["coordstring"]
Then that string can be used to initialize an astropy.coordinates.SkyCoord instance for each source.
Warning: The following cell takes a couple of minutes to run.
sdata["coord"] = coord.SkyCoord(sdata["coordstring"], unit=("hourangle", "deg"))
sdata["coord"]
And that in turn can be used to find the separation of each source from $\theta^1$ C.
sdata['Dprime'] = [c0.separation(c).arcsec for c in sdata["coord"]]
len(sdata['Dprime'])
rmax = 1200.0
nbins = rmax//4
_ = plt.hist(sdata["Dprime"], bins=nbins, range=(0, rmax), cumulative=True, color='b', alpha=0.3, label='Robberto 2013')
_ = plt.hist(sdata["Dprime"][m_rob_ionized], bins=nbins, range=(0, rmax), cumulative=True, color='r', alpha=0.3, label='Ionized, Robberto 2013')
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=True, color='y', alpha=0.7, label='Proplyds, Ricci 2008')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of stars with projected distance < D\'')
_ = plt.title('Cumulative distribution of stars vs distance: Robberto ACS sample')
plt.legend(loc='upper left')
plt.gcf().set_size_inches(18, 8)
rmax = 1200.0
binsize = 25
nbins = rmax//binsize
_ = plt.hist(sdata["Dprime"], bins=nbins, range=(0, rmax), cumulative=False, color='b', alpha=0.3, label='Robberto 2013')
_ = plt.hist(sdata["Dprime"][m_rob_ionized], bins=nbins, range=(0, rmax), cumulative=False, color='r', alpha=0.3, label='Ionized, Robberto 2013')
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=False, color='y', alpha=0.7, label='Proplyds, Ricci 2008')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of stars with projected distance D\' (25 arcsec bins)')
_ = plt.title('Distribution of stars vs distance: Robberto ACS sample')
plt.legend(loc='upper right')
plt.gcf().set_size_inches(18, 8)
rmax = 1200.0
sns.distplot(mirrored(sdata['Dprime']),
rug=True,
fit=stats.semicircular,
kde_kws={"kernel": "gau", "label": "KDE"},
rug_kws={"alpha": 0.03},
fit_kws={"label": "Semi-circular fit"},
bins=2*rmax//25,
label="Histogram (25 arcsec bins)", axlabel=r"Projected distance, $D'$")
plt.legend()
plt.xlim(0.0, 1200.0)
plt.title("Distribution of projected distances for all stars from Robberto et al. (2013) sample")
None
Zoom in on inner 200''
rmax = 200.0
nbins = rmax
_ = plt.hist(sdata["Dprime"], bins=nbins, range=(0, rmax), cumulative=True, color='b', alpha=0.3, label='Robberto 2013')
_ = plt.hist(sdata["Dprime"][m_rob_ionized], bins=nbins, range=(0, rmax), cumulative=True, color='r', alpha=0.3, label='Ionized, Robberto 2013')
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=True, color='y', alpha=0.7, label='Proplyds, Ricci 2008')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of stars with projected distance < D\'')
_ = plt.title('Cumulative distribution of stars vs distance: Robberto ACS sample')
plt.legend(loc='upper left')
plt.gcf().set_size_inches(18, 8)
Zoom in on inner 30''
rmax = 30.0
nbins = rmax
_ = plt.hist(sdata["Dprime"], bins=nbins, range=(0, rmax), cumulative=True, color='b', alpha=0.3, label='Robberto 2013')
_ = plt.hist(sdata["Dprime"][m_rob_ionized], bins=nbins, range=(0, rmax), cumulative=True, color='r', alpha=0.3, label='Ionized, Robberto 2013')
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=True, color='y', alpha=0.7, label='Proplyds, Ricci 2008')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of stars with projected distance < D\'')
_ = plt.title('Cumulative distribution of stars vs distance: Robberto ACS sample')
plt.legend(loc='upper left')
plt.gcf().set_size_inches(18, 8)
Maps of different source samples¶
rob_ra = [c.ra.deg for c in sdata['coord']]
rob_dec = [c.dec.deg for c in sdata['coord']]
rob_ra = np.array(rob_ra)
rob_dec = np.array(rob_dec)
ricci_ra = np.array([c.ra.deg for c in rdata['coord']])
ricci_dec = np.array([c.dec.deg for c in rdata['coord']])
import aplpy
import os
from imp import reload
aplpy.__version__
hstdir = os.path.expanduser("~/Dropbox/JorgeBowshocks/HST")
fig = aplpy.FITSFigure(os.path.join(hstdir, "mosaicf656-fixw-align.fits"))
fig.show_grayscale(vmin=0.0, vmax=5.e3)
fig._figure.set_size_inches(18, 18)
fig.show_markers(ricci_ra[m], ricci_dec[m], layer='Ricci proplyds', edgecolor='c', lw=3, alpha=0.8)
fig.show_markers(rob_ra[m_rob_ionized], rob_dec[m_rob_ionized], layer='Robb proplyds', edgecolor='orange')
fig.show_markers(rob_ra[~m_rob_ionized], rob_dec[~m_rob_ionized], layer='Robb non-proplyds', edgecolor='y', marker='*')
fig.recenter(c0.ra.deg, c0.dec.deg, radius=60/3600)
ax = fig._figure.axes[0]
ax.set_title('A title')
fig._figure
So the conclusion of this section is that the Ricci proplyd catalog (shown by blue circles on the image above) is more reliable than the Robberto ionized catalog (shown by orange circles) for the close-in proplyds. Robberto is missing a lot of obvious proplyds, and furthermore has a few sources that don't seem to correspond to anything, although some might be bowshocks or jet knots.
The Ricci catalog is only missing a very few, such as the one very close to th1c.
Monte Carlo simulations¶
Assume power law distribution in radius.
class SourcePopulation(object):
def __init__(self, N, Dmin=0.0, Dmax=800.0, m=2.6):
"""Return `N` samples of sources drawn from a radial distribution with density ~ D^-m
Maximum and minimum radii can be specified with `Dmin` and `Dmax`.
If `Dmin = 0` then `m < 3` is required."""
assert(m < 3)
self.N = N
self.m = m
self.Dmin = Dmin
self.Dmax = Dmax
# Radial distances from center
xi = np.random.random_sample((N,))
x0 = Dmin/Dmax
self.D = Dmax*((1 - x0**(3.0 - m))*xi + x0**(3.0 - m))**(1.0/(3.0 - m))
# Angle cosine with LOS
xi = np.random.random_sample((N,))
self.mu = 2*xi - 1.0
# Position angle
xi = np.random.random_sample((N,))
self.PA = 360.0*xi
# Inclination to plane of sky inc = [-90:90]
self.inc = np.degrees(np.arcsin(self.mu))
# Projected distance
self.Dprime = self.D*np.sqrt(1.0 - self.mu**2)
# Plane-of-sky coords
self.Dec = self.Dprime*np.cos(np.radians(self.PA))
self.RA = self.Dprime*np.sin(np.radians(self.PA))
p = SourcePopulation(150, m=2.1, Dmin=5.0, Dmax=250.0)
p2 = SourcePopulation(30, m=2.5, Dmin=150.0, Dmax=800.0)
Dprimes = np.r_[p.Dprime, p2.Dprime]
rmax = 800.0
nbins = rmax//4
_ = plt.hist(rdata["Dprime"][m], bins=nbins, range=(0, rmax), cumulative=True, color='b', alpha=0.3, label='Ricci 2008')
_ = plt.hist(Dprimes, bins=nbins, range=(0, rmax), cumulative=True, color='r', alpha=0.3, label='Fake sample')
# _ = plt.hist(p.Dprime, bins=nbins, range=(0, rmax), cumulative=True, color='g', alpha=0.3, label='p')
# _ = plt.hist(p2.Dprime, bins=nbins, range=(0, rmax), cumulative=True, color='y', alpha=0.8, label='p2')
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of proplyds with projected distance < D\'')
_ = plt.title('Cumulative distribution of proplyds vs distance: Both samples')
plt.legend(loc='upper left')
plt.gcf().set_size_inches(18, 8)
import seaborn
cmap = seaborn.palettes.dark_palette('red', as_cmap=True)
rmax = 800.0
nbins = rmax//25
datasets = [rdata["Dprime"][m], Dprimes, p.Dprime, p2.Dprime]
colors = ['b', 'y', 'g', 'y']
labels = ['Ricci 2008', 'Fake sample', 'p', 'p2']
_ = plt.hist(datasets, bins=nbins, range=(0, rmax), cumulative=False, color=None, label=labels)
_ = plt.xlabel('Projected distance D\', arcsec')
_ = plt.ylabel('# of proplyds with projected distance = D\' (bin width 25 arcsec)')
_ = plt.title('Distribution of proplyds vs distance: Both samples')
plt.legend(loc='upper right')
plt.gcf().set_size_inches(18, 8)
Dmax = 110.0
Dmaxx = 800.0
mzoom = Dprimes < Dmax
dataset = Dprimes[mzoom]
sns.distplot(dataset,
rug=True,
fit=stats.uniform,
kde_kws={"kernel": "gau", "label": "Gaussian KDE"},
rug_kws={"alpha": 0.3},
fit_kws={"label": "Constant fit"},
hist_kws={"range": (0, Dmaxx)},
bins=Dmaxx//10,
label="Histogram (10 arcsec bins)", axlabel=r"Projected distance, $D'$")
plt.legend(loc="upper left")
plt.xlim(0.0, 100.0)
plt.title("Distribution of projected distances from Monta Carlo fake sample");
Compare the projected separation and the inclination for the inner 30 arcsec of the MC sample¶
Assume that the inner wind shock lies at about $D = 26''$. Only proplyds that are physically clser than this will have Type I bowshocks. These are shown by orange points in the figure. Physically more distant proplyds are shown by gray points. We also show contours of physical separation every $10''$ from $D = 10''$ to $D = 200''$.
color = ".5"
m30 = p.Dprime <= 30.0
Dmax = 26.0
mwind = p.D <= Dmax
g = sns.JointGrid(p.Dprime[m30], p.inc[m30], size=9, xlim=(0, 30), ylim=(-90, 90))
g.plot_marginals(sns.distplot, color=color, rug=True,
kde_kws=dict(label="KDE", bw="silverman", cut=0))
#g.plot_joint(sns.regplot, color=color, fit_reg=False, scatter_kws=dict(s=50))
g.ax_joint.scatter(p.Dprime[mwind], p.inc[mwind], marker='o', c="orange", s=50, label="D < {}".format(Dmax))
g.ax_joint.scatter(p.Dprime[~mwind], p.inc[~mwind], marker='o', c=color, s=50, label="D > {}".format(Dmax))
inc_pts = np.linspace(-90.0, 90.0)
Dprime_pts = Dmax*np.cos(np.radians(inc_pts))
g.ax_joint.plot(Dprime_pts, inc_pts, '-')
for D in np.arange(10, 200, 10):
g.ax_joint.plot(Dprime_pts*D/Dmax, inc_pts, '-', color='k', lw=0.5, alpha=0.3)
g.ax_joint.legend(loc="upper center", frameon=True, shadow=True)
g.ax_joint.set_xlabel("Projected separation, $D'$")
g.ax_joint.set_ylabel("Inclination, degrees")
g.ax_marg_x.legend(loc="upper right")
g.ax_marg_y.legend(loc="best")
g.fig.tight_layout()
#g.fig.savefig("proplyd-joint-separation-size.pdf");
Distribution of sizes¶
Read in the size table that I converted from the DS9 regions.
sizedata = Table.read("Proplyd-Datasets/proplyd-chords-D60.dat", format='ascii.commented_header')
sizedata[:4]
Create coord.SkyCoord objects for the two ends of the chord that I measured at the back of each proplyd's head. Then use those to find the chord diameter in arcsec (Chord column), and then convert to a radius in units of $10^{14}$ cm (r14 column).
coords1 = coord.SkyCoord(ra=sizedata["RA1"], dec=sizedata["DEC1"], unit=("hourangle", "deg"))
coords2 = coord.SkyCoord(ra=sizedata["RA2"], dec=sizedata["DEC2"], unit=("hourangle", "deg"))
sizedata["Chord"] = [c1.separation(c2).arcsec for c1, c2 in zip(coords1, coords2)]
sizedata["r14"] = 0.5 * sizedata["Chord"] * 430*u.AU.in_units(u.cm)/1.e14
sizedata[:4]
Take the average of the RA and DEC components of the two ends of the chord to get an average coordinate (coords0) for each source. Note: This will not quite be the center of the proplyd, but it is close enough.
coords0 = coord.SkyCoord(ra=0.5*(coords1.ra + coords2.ra), dec=0.5*(coords1.dec + coords2.dec))
coords0[:4]
Remove the columns with the chord end coordinates because we don't need them any more.
sizedata.remove_columns(["RA1", "RA2", "DEC1", "DEC2"])
sizedata[:4]
Add some more columns:
RAandDeccorresponding tocoords0- Projected distance
Dprimefrom th1C (in arcsec) - Position angle
PAfrom th1C (in degrees).
sizedata['RA'] = coords0.ra.to_string(unit=u.hourangle, sep=':')
sizedata['Dec'] = coords0.dec.to_string(sep=':')
sizedata['Dprime'] = c0.separation(coords0).arcsec
sizedata['PA'] = c0.position_angle(coords0).deg
sizedata[-4:]
Now we try and merge in the HA98 table.
from astropy.table import join, Column
First naive attempt was a straight outer join on the ID column. This gave some objects that are in the HA98 table but not in the new data. Of these, 244-440 is OK, since it is outside our distance cut-off, but the others need fixing. 177-341 (not in Ricci catalog) is just a matter of the name, which should be 177-341W. There were also two really missing: 171-334 and 154-324 (my typo).
So we need to fix things up. First make the ID field in the HA98 table wider and fix the name of 177-341W.
This took much longer than it should have done. It turns out that you can't change the dtype of a column without making a new table.
data = Table(data, names=data.colnames, dtype=['<U12', '<f8', '<f8', '<f8'])
data['ID'] = Column(data['ID'], dtype='<U12')
irow = np.argmax(data['ID'] == '177-341')
data['ID'][irow] = '177-341W'
So now we can do the join properly:
jdata = join(sizedata, data, join_type='outer', keys='ID')
And the only one missing is 244-440:
jdata[~jdata["r14_2"].mask & jdata["r14_1"].mask]
Now we tidy up some of the column names:
jdata.rename_column("r14_1", "r14")
jdata.rename_column("r14_2", "r14_HA")
jdata.rename_column("d", "Dprime_HA")
Look at the rows that are in both datasets:
overlap = ~jdata["r14"].mask & ~jdata["r14_HA"].mask
jdata[overlap]
Check that the two estimates of the proplyd sizes are not too different.
x, y = [np.log10(jdata[overlap][colname]) for colname in ['r14', 'r14_HA']]
g = sns.jointplot(x, y, size=8)
g.ax_joint.plot([0, 2], [0, 2], '-w', zorder=-100, lw=1)
g.ax_joint.set_xlim(0.0, 2.0)
g.ax_joint.set_ylim(0.0, 2.0);
So it doesn't look too bad. Here are the cases where HA98 got >30% larger sizes. These are mainly large sized.
jdata[overlap][jdata[overlap]['r14_HA'] > 1.3*jdata[overlap]['r14']]['ID', 'r14', 'r14_HA']
And here are the cases where HA98 sizes were more than 30% smaller. These are mainly small sized. In this case, I think that the current results are more reliable because in HA98 we did not remove the central stars and so overestimated the compactness of some of the small proplyds.
jdata[overlap][jdata[overlap]['r14'] > 1.3*jdata[overlap]['r14_HA']]['ID', 'r14', 'r14_HA']
If we look at the distribution of ratios, then it seems we have agreement to a level of about 30%. This is not bad, given that the sizes range over an order of magnitude.
ratio = jdata['r14']/jdata['r14_HA']
sns.distplot(ratio[overlap], rug=True, fit=stats.norm, axlabel='r14 / r14_HA')
plt.title('Current vs HA98 estimates of proplyd size:'
' ratio = {:.3f} +/- {:.3f}'.format(ratio[overlap].mean(),
ratio[overlap].std()))
We can do the same thing for the separations. This shows a much tighter correlation, as one would expect
x, y = [np.log10(jdata[overlap][colname]) for colname in ['Dprime', 'Dprime_HA']]
g = sns.jointplot(x, y, size=8)
g.ax_joint.plot([0, 2], [0, 2], '-w', zorder=-100, lw=1)
g.ax_joint.set_xlim(0.0, 2.0)
g.ax_joint.set_ylim(0.0, 2.0);
And just for fun, here is the same plot for the Vicente & Alves (2005) sizes.
x, y = [np.log10(jdata[overlap][colname]) for colname in ['VicR14', 'r14']]
g = sns.jointplot(x, y, color='brown', size=8)
g.ax_joint.plot([0, 2], [0, 2], '-w', zorder=-100, lw=1)
g.ax_joint.set_xlim(0.0, 2.0)
g.ax_joint.set_ylim(0.0, 2.0);
Finally, we can go back to the size-distance relationship. First, the joint distribution as a scatter plot, with the marginal distributions shown on the top and right edges.
color = "#774499"
# g = sns.JointGrid(mirrored(jdata["Dprime"].filled(0.0)), doubled(jdata["r14"].filled(0.0)),
g = sns.JointGrid(jdata["Dprime"], jdata["r14"],
size=9, xlim=(0, 60), ylim=(0, 35))
# g.plot(sns.regplot, sns.distplot, stats.pearsonr);
g.plot_marginals(sns.distplot, color=color, bins=10, rug=True, kde=False)
g.plot_joint(sns.regplot, color=color, ci=99, order=2, label="Quadratic regression fit with 99% confidence bounds")
g.ax_joint.legend(loc="upper left")
g.ax_joint.set_xlabel("Projected separation, $D'$")
g.ax_joint.set_ylabel("Proplyd radius $r_0$, $10^{14}$ cm")
g.fig.tight_layout()
g.fig.savefig("proplyd-joint-separation-size-new.pdf");
Now, split up into coarse bins according to the separation, and then look at the distribution of sizes in each bin.
First, sort on separation (it turns out that sorting doesn't work for colums that are masked arrays, so we fill the missing values with NaNs first):
jjdata = jdata[:-1].filled(np.nan)
jjdata.sort('Dprime')
Now split into 5 equal parts and make a violin plot of each one:
nsplit = 5
Dprime_splits = np.array_split(jjdata['Dprime'], nsplit)
r14_splits = np.array_split(jjdata['r14']*1e14*u.cm.to(u.au), nsplit)
split_names = ['D\' = {:.1f} to {:.1f}"'.format(x.min(), x.max()) for x in Dprime_splits]
sns.violinplot(r14_splits,
names=split_names,
inner="points", color="PuBuGn")
plt.ylim(0.0, None)
plt.ylabel("Proplyd size, AU")
plt.xlabel("Range of projected separations");
There are 14 (or 13) sources in each bin.
[len(a) for a in r14_splits]
So it looks like there is no significant change in the small proplyds with $r < 10^{15}$ cm (60 AU). The larger proplyds are reduced in number in the inner 3 bins.
So, for a first stab at the size distribution, I will try and fit something sensible to the 27 sources with $D' > 35.5$. And then play with having a varying upper cut-off.
But first, we add a new column with radius in AU, and then save the whole table to a file so that we can use it later more easily.
jjdata['rAU'] = jjdata['r14']*(1.e14*u.cm).to(u.AU)
jjdata[:4]
jjdata.write("Proplyd-datasets/merged-proplyd-sizes.tab",
format="ascii.tab",
fill_values=[('nan', '--')])
outer = jjdata['Dprime'] > 30.0
plt.hist([jjdata['rAU'], jjdata[~outer]['rAU']],
label=['All', 'D\' > 35"'],
normed=True, cumulative=-1, bins=50)
plt.legend()
x = np.linspace(0.0, jjdata['rAU'].max())
r0, a = 50.0, 1.1
rmax = 150.0
b = 1.0 - 1.0/(1.0 + (rmax/r0)**(2*a))
plt.plot(x, (1.0/(1.0 + (x/r0)**(2*a)) - (1.0 - b))/b, 'w')
# plt.plot(x, , 'w')
plt.xlabel('Proplyd radius, r, AU')
plt.ylabel('Fraction of sources smaller or equal to r')
None
So it looks like a reasonable fit to the CDF is $1 - 1/(1 + (r/r_0)^{2a})$ with $r_0 = 55''$ and $a = 1.3$.
Which means that the PDF is $\sim r^{2a - 1} / (1 + (r/r_0)^{2a})^2$.
nbins = 15
plt.hist([jjdata['rAU'], doubled(jjdata[outer]['rAU'])],
label=['All', 'D\' > 35"'],
normed=False, bins=nbins)
plt.legend()
x = np.linspace(0.0, jjdata['rAU'].max())
r0, a = 70.0, 1.1
plt.plot(x, 5*(len(jjdata)/nbins)*2*a*(x/r0)**(2*a - 1)/(1.0 + (x/r0)**(2*a))**2, 'w')
# plt.plot(x, , 'w')
plt.xlabel('Proplyd radius, r, AU')
plt.ylabel('PDF');
So that looks sort of alright, although it might be said to peak at too large a value.
So we can do the random sampling from a uniform deviate $\xi$, where $r = r_0 (\xi^{-1} - 1)^{1/(2a)}$. That is assuming no cut-off. I changed the parameters in the end to $r_0 = 50''$ and $a = 1.4$.
If there is a cut-off for $r > r_1$, then we can define $b = 1 - \left(1 + (r_1/r_0)^{2a}\right)^{-1}$ and then $r = r_0 \left((1 - b \xi)^{-1} - 1\right)^{1/(2a)}$.
r0, a = 70.0, 1.2
rmax = 300.0
b = 1.0 - 1.0/(1.0 + (rmax/r0)**(2*a))
print(b)
samples = [jjdata['rAU']]
names = ["proplyds"]
for i in range(4):
xi = np.random.random_sample((len(jjdata),))
samples.append(r0*((1.0/(1.0 - b*xi)) - 1.0)**(1./(2*a)))
names.append('fake #'+str(i+1))
sns.violinplot(samples, names=names,
inner="points", color="PuBuGn")
plt.ylim(0.0, 350.0)
plt.title("Proplyd size distribution: real versus simulated")
plt.ylabel("Proplyd size, AU")
plt.xlabel("Sample");
That looks great. So now we can start worrying about the high-end cut-off and how it varies with separation.
r0, a = 50.0, 1.1
samples = []
names = []
for rmax in [100.0, 150.0, 250.0, 350.0, 350.0]:
b = 1.0 - 1.0/(1.0 + (rmax/r0)**(2*a))
xi = np.random.random_sample((14,))
samples.append(r0*((1.0/(1.0 - b*xi)) - 1.0)**(1./(2*a)))
names.append('fake, rmax = {:.0f}'.format(rmax))
sns.violinplot(samples, names=names,
inner="points", color="PuBuGn")
plt.ylim(0.0, 350.0)
plt.title("Proplyd size distribution: simulated with different cut-offs rmax")
plt.ylabel("Proplyd size, AU")
plt.xlabel("Sample");
So now it looks like a shallower power law is necessary: $r_0 = 50''$, $a = 1.1$. Otherwise, there are not enough sources in the high tail of the size distributions and the cut-off doesn't get a chance to do anything.
Try adding these samples together, to see how well it fits the observed total distro.
total_fake = np.hstack(samples)
sns.violinplot([jjdata['rAU'], total_fake], names=['proplyds', 'simulated'],
inner="points", color="PuBuGn")
plt.ylim(0.0, 350.0)
plt.title("Proplyd size distribution: real vs simulated using 4 different rmax")
plt.ylabel("Proplyd size, AU")
plt.xlabel("Sample");
So that is sort of OK. But we really need to have the upper cut-off be a function of true separation $D$. Then we need to choose the $D$'s from a distribution, like we did earlier.
First of all, decide what is the functional form of $r_1$ versus $D$. We can check the values of $D$ for sources in our previous random sample p that have $D'$ values that lie in each of our Dprime_splits.
D_pts = []
for Dp in Dprime_splits:
mask = (p.Dprime > Dp.min()) & (p.Dprime <= Dp.max())
DD = p.D[mask]
print(DD.min(), np.median(DD), DD.max())
D_pts.append(np.median(DD))
rmax_pts = [100, 150, 250, 350, 350]
plt.plot(D_pts, rmax_pts, 'o')
D_ppts = np.r_[np.linspace(0.0, 0.99, 100), np.linspace(1.0, 100.0, 200)]
plt.plot(D_ppts, 40*np.sqrt(D_ppts))
plt.plot(D_ppts, 6*D_ppts)
plt.ylim(0.0, 400.0)
plt.xlabel('True separation, $D$, arcsec')
plt.ylabel('Upper cut-off on proplyd size, $r_1$, AU');
So it looks like $r_1 = 40 (D\,/\,'')^{1/2}$ AU would work well. So we will now try and put it all together.
class SourcePopulation(object):
def __init__(self, N, Dmin=0.0, Dmax=800.0, m=2.6, r0=50.0, a=1.1):
"""Return `N` samples of proplyds drawn from a radial distribution with density ~ D^-m
Maximum and minimum radii can be specified with `Dmin` and `Dmax`.
If `Dmin = 0` then `m < 3` is required. The proplyd size distribution
is described by `r0` and `a`"""
assert(m < 3)
self.N = N
self.m = m
self.Dmin = Dmin
self.Dmax = Dmax
# Radial distances from center
xi = np.random.random_sample((N,))
x0 = Dmin/Dmax
self.D = Dmax*((1 - x0**(3.0 - m))*xi + x0**(3.0 - m))**(1.0/(3.0 - m))
# Angle cosine with LOS
xi = np.random.random_sample((N,))
self.mu = 2*xi - 1.0
# Position angle
xi = np.random.random_sample((N,))
self.PA = 360.0*xi
# Inclination to plane of sky inc = [-90:90]
self.inc = np.degrees(np.arcsin(self.mu))
# Projected distance
self.Dprime = self.D*np.sqrt(1.0 - self.mu**2)
# Plane-of-sky coords
self.Dec = self.Dprime*np.cos(np.radians(self.PA))
self.RA = self.Dprime*np.sin(np.radians(self.PA))
# Proplyd size in AU
# rmax = 40*np.sqrt(self.D)
rmax = 6*self.D
b = 1.0 - 1.0/(1.0 + (rmax/r0)**(2*a))
xi = np.random.random_sample((N,))
self.rAU = r0*((1.0/(1.0 - b*xi)) - 1.0)**(1./(2*a))
num60 = []
mega_datasets = []
for _ in range(30000):
p = SourcePopulation(190, m=2.1, Dmin=5.0, Dmax=250.0)
num60.append(np.sum(p.Dprime < 60.0))
mega_datasets.append(p)
We play with the total number of proplyds until we get the right number of proplyds at projected distances less than $60''$, which should be 69. Looks like 190 is about right.
sns.distplot(num60, bins=20, kde=False, fit=stats.norm);
#istart = 1133
#istart = 7777
istart = 999
scatter_kws = dict(marker='o', alpha=1.0, edgecolors=None)
line_kws = dict(lw=2)
fig, axes = plt.subplots(5, 5, sharex=True, sharey=True)
proplyd_colors=sns.color_palette("RdPu", 5)
sns.regplot(jjdata['Dprime'], jjdata['rAU'], ax=axes[0, 0], order=2, ci=99,
color=proplyd_colors[3], scatter_kws=scatter_kws, line_kws=line_kws)
axes[0, 0].text(3, 340, 'Orion proplyds, N = '+str(len(jjdata)))
axes[0, 0].set_xlabel('')
i = istart
for p, ax in zip(mega_datasets[istart:istart+24], axes.flat[1:]):
m60 = p.Dprime <= 60.0
color_kws = scatter_kws.copy()
color_kws.update(c=np.abs(p.inc[m60]), cmap="PuBuGn", vmin=0.0, vmax=90.0)
g = sns.regplot(p.Dprime[m60], p.rAU[m60], ax=ax, order=2, ci=99,
color="#774499", scatter_kws=color_kws, line_kws=line_kws)
ax.text(3, 340, 'Simulation #{}, N = {}'.format(i, m60.sum()))
i += 1
for ax in axes[-1, :]:
ax.set_xlabel("$D'$, arcsec")
for ax in axes[:, 0]:
ax.set_ylabel('$r$, AU')
for ax in axes.flat:
ax.tick_params(axis='both', which='major', labelsize=9)
ax.set_xlim(0.0, 62.0)
ax.set_ylim(0.0, 370.0)
cb = plt.colorbar(g.collections[0], ax=list(axes.flat), shrink=0.25, fraction=0.02)
cb.set_label("Inclination, deg")
fig.suptitle("Comparison of Orion proplyds with 24 simulated populations", fontsize='xx-large')
#fig.tight_layout(rect=(0, 0, 1, 0.98), h_pad=0.1, w_pad=0.1)
fig.subplots_adjust(left=0.1, right=0.87, bottom=0.1, top=0.96, wspace=0.05, hspace=0.05)
fig.set_size_inches(16, 16)
These work out great. The majority of the simulated datasets bear at least a passing resemblance to te real one. I had to change the cut-off size to be linear in $D$, otherwise there were too many simulated datasets with large proplyds at small distances.
For the simulated populations we show the absolute value of the inclination. Obviously, we can't know this for the proplyds (actually, maybe we could).
Note that the large proplyds at small radii tend to be at high inclinations, meaning that the true $D$ is much larger than the projected $D'$.
#istart = 1133
#istart = 7777
istart = 999
nsplit = 5
scatter_kws = dict(marker='o', alpha=0.6, edgecolors='none')
line_kws = dict(lw=2)
fig, axes = plt.subplots(5, 5, sharex=False, sharey=True)
# Plot real proplyd dataset
rAU_splits = np.array_split(jjdata['rAU'], nsplit)
Dp_splits = np.array_split(jjdata['Dprime'], nsplit)
labels = ['{:.0f}–{:.0f}'.format(x.min(), x.max()) for x in Dp_splits]
sns.violinplot(rAU_splits, inner="points", names=labels, color="RdPu", ax=axes[0, 0])
axes[0,0].set_xlabel('')
# Plot a bunch of fake datasets
for p, ax in zip(mega_datasets[istart:istart+24], axes.flat[1:]):
m60 = p.Dprime <= 60.0
sort_order = np.argsort(p.Dprime[m60])
rAU_splits = np.array_split(p.rAU[m60][sort_order], nsplit)
Dp_splits = np.array_split(p.Dprime[m60][sort_order], nsplit)
labels = ['{:.0f}–{:.0f}'.format(x.min(), x.max()) for x in Dp_splits]
sns.violinplot(rAU_splits, inner="points", names=labels, color="PuBuGn", ax=ax)
# Sort out the axis labels
for ax in axes[:, 0]:
ax.set_ylabel('$r$, AU')
for ax in axes[-1, :]:
ax.set_xlabel("Ranges of $D'$, arcsec")
for ax in axes.flat:
ax.tick_params(axis='both', which='major', labelsize=9)
ax.set_ylim(0.0, 370.0)
fig.suptitle("Comparison of Orion proplyds with 24 simulated populations", fontsize='xx-large')
fig.tight_layout(rect=(0, 0, 1, 0.97), h_pad=0.2, w_pad=0.2)
fig.set_size_inches(16, 16)
The distribution of $\beta$¶
Finally we can get on to this. First we work out the density at the i-front on the axis and thus the equivalent spherical $\dot{M}$. Note that this is larger than the true mass-loss rate, since it doesn't account for the reduction in density away from the axis.
from astropy.constants import c, m_p, Ryd, L_sun, M_sun
from astropy.constants import h as hplanck
QH = 1e49/u.s
d_orion = 430
alpha_B = 2.6e-13*u.cm**3/u.s
alpha_ha = 1.0e-13*u.cm**3/u.s
omega = 0.1
# ionized sound speed
ci = (11.0*u.km/u.s).to(u.cm/u.s)
# mean mass per particle
m = (1.3*u.u).to(u.g)
# H alpha photon energy
E_ha = (hplanck*c/(6563.0*u.angstrom)).to(u.erg)
#iset = 1021
iset = 1999
p = mega_datasets[iset]
D = (p.D*d_orion*u.au).to(u.cm)
F = QH/(4.0*np.pi*D**2)
r = (p.rAU*u.AU).to(u.cm)
# u0 n + alpha_B omega n^2 r = F
ntilde = ci/(2*alpha_B*omega*r)
n = np.sqrt(ntilde**2 + F / (alpha_B*omega*r)) - ntilde
lambda_ = ci*n/F
L_ha = E_ha*alpha_ha*r**3*n**2
S_ha = E_ha*alpha_ha*r*n**2/(4.0*np.pi*u.steradian)
#lambda_[(p.Dprime < 60.0) & (p.rAU > 50.0)].max()
(L_ha/L_sun).decompose().max()
Mdot = (n*ci*m*4*np.pi*r**2).to(u.solMass/u.year)
Mdot.mean()
We look at the distribution of the advection parameter $\lambda = n_0 u_0 / F$ versus the H$\alpha$ luminosity.
m60 = p.Dprime < 60.0
x = L_ha[m60]
y = lambda_[m60]
g = sns.jointplot(np.log10(x.value), np.log10(y))
g.ax_joint.set_ylabel(r'Advection parameter: $\log_{10} \,\lambda$')
#g.ax_joint.set_xlabel('Proplyd luminosity $\\log_{10}L(\\mathrm{H\\alpha})$ / [{:LATEX}])'.format(x.unit));
g.ax_joint.set_xlabel('Proplyd luminosity: $\\log_{{10}}\\,L$(H$\\alpha$) [{:LATEX}]'.format(x.unit))
And the surface brightness versus distance:
rootcolor = (0.5, 0.7, 0.9)
ncolors = 16
darkmap = sns.dark_palette(rootcolor, ncolors, reverse=True, as_cmap=True)
sns.palplot(sns.dark_palette(rootcolor, ncolors, reverse=True))
First we use a jointplot to see the marginal distributions.
x = p.Dprime[m60]
y = S_ha[m60]
z = p.rAU[m60]
inc = np.abs(p.inc[m60])
g = sns.jointplot(np.log10(x), np.log10(y.value),
joint_kws={"s": z, "c": inc, "vmin": 0.0, "vmax": 90.0, "cmap": darkmap},
size=8)
Then we try with regplot, plotting a linear regression. The results vary a bit between samples, but generally the regression line is a bit steeper than $D'^{-2}$ because the average inclination becomes higher at smaller separations.
This compares pretty well with O'Dell (1998) Fig 2, once we take into account that he is using photon units.
g = sns.regplot(np.log10(x), np.log10(y.value), order=1,
robust=False,
scatter_kws={"s": z, "c": inc, "vmin": 0.0, "vmax": 90.0, "cmap": darkmap}
)
g.set_ylim(-1.0, 3.0)
xpts = np.linspace(*g.get_xlim())
for y0 in [2.0, 3.0, 4.0, 5.0]:
g.plot(xpts, y0 - 2*xpts, '-', c='w', zorder=-100, lw=1)
cb = plt.colorbar(g.collections[0], ax=g)
cb.set_label("Inclination, deg")
g.set_xlabel("log10(D' / arcsec)")
g.set_ylabel("H alpha brightness, log10(S / erg/s/cm^2/sr)")
g.set_title("Surface brightness versus projected separation")
plt.gcf().set_size_inches(10, 8)
We can get a better regression by controlling for inclination as a confounding variable.
g = sns.regplot(np.log10(x), np.log10(y.value), order=3,
x_partial=1.0/np.cos(np.radians(inc)),
y_partial=1.0/np.cos(np.radians(inc)),
# robust=True,
scatter_kws={"s": z, "c": inc, "vmin": 0.0, "vmax": 90.0, "cmap": darkmap}
)
g.set_ylim(-1.0, 3.0)
xpts = np.linspace(*g.get_xlim())
for y0 in [2.0, 3.0, 4.0, 5.0]:
g.plot(xpts, y0 - 2*xpts, '-', c='w', zorder=-100, lw=1)
cb = plt.colorbar(g.collections[0], ax=g)
cb.set_label("Inclination, deg")
g.set_xlabel("log10(D' / arcsec)")
g.set_ylabel("H alpha brightness, log10(S / erg/s/cm^2/sr)")
g.set_title("Controlling for inclination as a confounding variable")
plt.gcf().set_size_inches(10, 8)
#m60 = np.ones_like(p.Dprime).astype(bool)
m60 = p.Dprime < 60
x = p.Dprime[m60]
y = Mdot[m60]
z = p.rAU[m60]
inc = np.abs(p.inc[m60])
g = sns.jointplot(np.log10(x), np.log10(y.value),
joint_kws={"s": z, "c": inc, "vmin": 0.0, "vmax": 90.0, "cmap": darkmap},
size=8)
Here is the mass loss rate distribution for a bunch of simulated populations. There is a general tendency for the regression line to fall slightly with projected distance. But for the highest mass-loss rates, there doesn't seem to be much of a trend with distance at all.
#istart = 999
istart = 1999
scatter_kws = dict(marker='o', alpha=1.0, edgecolors=None)
line_kws = dict(lw=2)
fig, axes = plt.subplots(5, 5, sharex=True, sharey=True)
i = istart
for p, ax in zip(mega_datasets[istart:istart+25], axes.flat):
m60 = p.Dprime <= 60.0
D = (p.D*d_orion*u.au).to(u.cm)
F = QH/(4.0*np.pi*D**2)
r = (p.rAU*u.AU).to(u.cm)
ntilde = ci/(2*alpha_B*omega*r)
n = np.sqrt(ntilde**2 + F / (alpha_B*omega*r)) - ntilde
Mdot = (n*ci*m*4*np.pi*r**2).to(u.solMass/u.year)
color_kws = scatter_kws.copy()
color_kws.update(c=np.abs(p.inc[m60]), s=0.5*p.rAU[m60], cmap="PuBuGn", vmin=0.0, vmax=90.0)
g = sns.regplot(p.Dprime[m60], Mdot[m60]/1.e-7, ax=ax, order=2, ci=99,
color="#774499", scatter_kws=color_kws, line_kws=line_kws)
ax.text(0.02, 0.92, 'Simulation #{}, N = {}'.format(i, m60.sum()), transform=ax.transAxes)
i += 1
for ax in axes[-1, :]:
ax.set_xlabel("$D'$, arcsec")
for ax in axes[:, 0]:
ax.set_ylabel(r'$\dot{M}$, 1e-7 solar masses / yr')
for ax in axes.flat:
ax.tick_params(axis='both', which='major', labelsize=9)
ax.set_xlim(0.0, 62.0)
ax.set_ylim(-0.3, 20.0)
cb = plt.colorbar(g.collections[0], ax=list(axes.flat), shrink=0.25, fraction=0.02)
cb.set_label("Inclination, deg")
fig.suptitle("Mass loss rates for 25 simulated populations", fontsize='xx-large')
#fig.tight_layout(rect=(0, 0, 1, 0.98), h_pad=0.1, w_pad=0.1)
fig.subplots_adjust(left=0.1, right=0.87, bottom=0.1, top=0.96, wspace=0.05, hspace=0.05)
fig.set_size_inches(16, 16)
Next, we do the same for the $\beta$ values, which are directly proportional so the graphs look pretty identical.
#istart = 999
istart = 1999
# proplyd flow Mach number at shock
Mach = 3.0
# O-star wind parameters
Mdw = 3.5e-7*u.solMass/u.year
Vw = (1200*u.km/u.s).to(u.cm/u.s)
scatter_kws = dict(marker='o', alpha=1.0, edgecolors=None)
line_kws = dict(lw=2)
fig, axes = plt.subplots(5, 5, sharex=True, sharey=True)
i = istart
for p, ax in zip(mega_datasets[istart:istart+25], axes.flat):
m60 = p.Dprime <= 60.0
D = (p.D*d_orion*u.au).to(u.cm)
F = QH/(4.0*np.pi*D**2)
r = (p.rAU*u.AU).to(u.cm)
ntilde = ci/(2*alpha_B*omega*r)
n = np.sqrt(ntilde**2 + F / (alpha_B*omega*r)) - ntilde
Mdot = (n*ci*m*4*np.pi*r**2).to(u.solMass/u.year)
beta = Mdot*Mach*ci/(Mdw*Vw)
color_kws = scatter_kws.copy()
color_kws.update(c=np.abs(p.inc[m60]), s=0.5*p.rAU[m60], cmap="PuBuGn", vmin=0.0, vmax=90.0)
g = sns.regplot(p.Dprime[m60], beta[m60], ax=ax, order=2, ci=99,
color="#774499", scatter_kws=color_kws, line_kws=line_kws)
ax.text(0.02, 0.92, 'Simulation #{}, N = {}'.format(i, m60.sum()), transform=ax.transAxes)
i += 1
for ax in axes[-1, :]:
ax.set_xlabel("$D'$, arcsec")
for ax in axes[:, 0]:
ax.set_ylabel(r'$\beta$')
for ax in axes.flat:
ax.tick_params(axis='both', which='major', labelsize=9)
ax.set_xlim(0.0, 62.0)
ax.set_ylim(-0.003, 0.15)
cb = plt.colorbar(g.collections[0], ax=list(axes.flat), shrink=0.25, fraction=0.02)
cb.set_label("Inclination, deg")
fig.suptitle("Wind-wind beta values for 25 simulated populations", fontsize='xx-large')
#fig.tight_layout(rect=(0, 0, 1, 0.98), h_pad=0.1, w_pad=0.1)
fig.subplots_adjust(left=0.1, right=0.87, bottom=0.1, top=0.96, wspace=0.05, hspace=0.05)
fig.set_size_inches(16, 16)
We can sum all the simulation samples to get better statistics
multiplier = 10000
bigp = SourcePopulation(190*multiplier, m=2.1, Dmin=5.0, Dmax=250.0)
Joint distribution of size and projected separation for the large sample. We can clearly see the corners due to the sharp cut-offs.
sns.jointplot(bigp.Dprime, bigp.rAU, stat_func=None, kind='hex', size=10, xlim=(0, 60), ylim=(0, 350))
m60 = bigp.Dprime <= 60.0
m30 = bigp.Dprime <= 30.0
mm30 = bigp.D <= 30.0
mm15 = bigp.D <= 15.0
D = (bigp.D*d_orion*u.au).to(u.cm)
F = QH/(4.0*np.pi*D**2)
r = (bigp.rAU*u.AU).to(u.cm)
ntilde = ci/(2*alpha_B*omega*r)
n = np.sqrt(ntilde**2 + F / (alpha_B*omega*r)) - ntilde
Mdot = (n*ci*m*4*np.pi*r**2).to(u.solMass/u.year)
beta = Mdot*Mach*ci/(Mdw*Vw)
Histograms of the log of the $\beta$ values for the different subsamples. The 60 arcsec sample tends to have slightly smaller values than the 30 arcsec sample. Also shown is the sample with true physical separation $< 30''$. This would be relevant if that was the physical extent of the hypersonic wind region. This last sample has a narrow distribution of $\beta$, missing both the high values and the low values. For instance, it drops to zero at about $\beta = 0.1$.
hist_kws = {"range": (-5.0, -0.5), "alpha": 0.6}
nbins = 100
#with sns.color_palette("RdPu_r", 4):
with sns.color_palette("coolwarm", 4):
ax = sns.distplot(np.log10(beta[m60]), kde=False, bins=nbins, hist_kws=hist_kws, label="D' <= 60")
sns.distplot(np.log10(beta[m30]), kde=False, bins=nbins, hist_kws=hist_kws, ax=ax, label="D' <= 30")
sns.distplot(np.log10(beta[mm30 & ~mm15]), kde=False, bins=nbins, hist_kws=hist_kws, ax=ax, label="15 < D <= 30")
sns.distplot(np.log10(beta[mm15]), kde=False, bins=nbins, hist_kws=hist_kws, ax=ax, label="D <= 15")
ax.set_xlim(*hist_kws["range"])
ax.legend()
ax.set_xlabel('log10(beta)')
ax.set_ylabel('PDF');
In the following two plots we show the joint distribution of inclination and $\beta$ for three different subsamples of the $D' <= 30''$ sample:
- True distance $D <= 15''$. The inner core, which is definitely within the supersonic wind region.. This has a rather flat distribution in $\beta$ that peaks at 0.01–0.04 and seems almost independent of inclination. The inclination distribition is roughly uniform in $\sin i$, as expected.
- True distance with $15'' < D <= 30''$. The intermediate shell, which may or may not be in the supersonic wind zone. These still have the same $i$ distribution as the inner core sample, but the $\beta$ distribution peaks at smaller values of 0.005–0.02.
- Projected distance $D' <= 30''$ but $D > 30''$. The outer interlopers, which should not show any wind-wind bowshocks since they are simply projected onto the inner core.
mask = mm15
g = sns.jointplot(np.abs(bigp.inc[mask]), beta[mask],
stat_func=None, kind='hex', size=6, xlim=(0, 90), ylim=(0, 0.1))
g.set_axis_labels('inclination', 'beta')
nexpect = mask.sum()/multiplier
g.ax_joint.text(2, 0.095,
"Sample with true separation $D < 15''$, $N = {:.1f}$".format(nexpect),
fontsize='x-large');
mask = mm30 & ~mm15
g = sns.jointplot(np.abs(bigp.inc[mask]), beta[mask],
stat_func=None, kind='hex', size=6, xlim=(0, 90), ylim=(0, 0.1))
g.set_axis_labels('inclination', 'beta')
nexpect = mask.sum()/multiplier
g.ax_joint.text(2, 0.095,
"Sample with true separation $15'' < D < 30''$, $N = {:.1f}$".format(nexpect),
fontsize='x-large');
mask = m30 & ~mm30
g = sns.jointplot(np.abs(bigp.inc[mask]), beta[mask],
stat_func=None, kind='hex', size=6, xlim=(0, 90), ylim=(0, 0.1))
g.set_axis_labels('inclination', 'beta')
nexpect = mask.sum()/multiplier
g.ax_joint.text(5, 0.09,
"Sample with $D > 30''$ but $D' < 30''$, $N = {:.1f}$".format(nexpect),
fontsize='x-large');
The next three graphs show the $\beta$ distribution in more detail for the three sub-samples. Separate histograms are shown for low and high inclinations. For the two inner sub-samples there is no difference between the two inclination ranges, which must be true since the selection criteria are independent of $i$. For the outer interlopers, on the other hand, there is a difference: smaller $\beta$ for higher inclinations, due to the implicit dependence of $D$ on $i$ for this sub-sample that is caused by imposing a cut-off in $D'$.
incmask = np.abs(bigp.inc) < 30.0
hist_kws = {"range": (-0.001, 0.101), "alpha": 0.4}
nbins = 100
with sns.color_palette("coolwarm", 2):
ax = sns.distplot(beta[mm15 & incmask], kde=False, bins=nbins, hist_kws=hist_kws, label="|i| < 30 deg")
sns.distplot(beta[mm15 & ~incmask], kde=False, bins=nbins, hist_kws=hist_kws, ax=ax, label="|i| >= 30 deg")
ax.set_xlim(*hist_kws["range"])
ax.legend()
ax.set_title('Inner core: D <= 15\'\'')
ax.set_xlabel('beta')
ax.set_ylabel('PDF');
with sns.color_palette("coolwarm", 2):
ax = sns.distplot(beta[mm30 & ~mm15 & incmask], kde=False, bins=nbins, hist_kws=hist_kws, label="|i| < 30 deg")
sns.distplot(beta[mm30 & ~mm15 & ~incmask], kde=False, bins=nbins, hist_kws=hist_kws, ax=ax, label="|i| >= 30 deg")
ax.set_xlim(*hist_kws["range"])
ax.legend()
ax.set_title('Intermediate shell: D = 15\'\'–30\'\'')
ax.set_xlabel('beta')
ax.set_ylabel('PDF');
incmask = np.abs(bigp.inc) < 60.0
with sns.color_palette("coolwarm", 2):
ax = sns.distplot(beta[m30 & ~mm30 & incmask], kde=False, bins=nbins, hist_kws=hist_kws, label="|i| < 60 deg")
sns.distplot(beta[m30 & ~mm30 & ~incmask], kde=False, bins=nbins, hist_kws=hist_kws, ax=ax, label="|i| >= 60 deg")
ax.set_xlim(*hist_kws["range"])
ax.legend()
ax.set_title('Outer interlopers: D\' < 30\'\' but D > 30\'\'')
ax.set_xlabel('beta')
ax.set_ylabel('PDF');
sns.distplot(np.sin(np.radians(bigp.inc[mm30])), kde=False, bins=nbins)
Finally, bowshock shapes¶
import sys
sys.path.insert(0, "./conic-projection")
import conproj_utils
conic = conproj_utils.Conic(A=2.0, th_conic=15.0)
trange = conic.make_t_array([conic.tparallel(0.0), np.pi/2])
plt.plot(conic.x(trange), conic.y(trange))
plt.axis('equal')
plt.xlim(-3.0, 1.2)
def A(beta):
return 1.5/(1.0 - np.sqrt(beta))
def thc(beta, xi=1.0):
arg = 3*(1.0/(1.0 - np.sqrt(beta)) - xi*(1.0 + np.sqrt(beta))**2/(1.0 - xi*beta)**2/(1 + 0.2*xi*beta))
return np.sign(arg)*np.arctan(np.sqrt(np.abs(arg)))
beta_grid = np.r_[np.linspace(0.0, 0.01, 100), np.linspace(0.011, 0.1, 100)]
with sns.color_palette("PuRd_d", 4):
plt.plot(beta_grid, A(beta_grid), label='A = R_c/R_0')
plt.plot(beta_grid, thc(beta_grid), label='theta_c (isotropic)')
plt.plot(beta_grid, thc(beta_grid, xi=0.8), label='theta_c (proplyd)')
plt.plot(beta_grid, thc(beta_grid, xi=0.6), label='theta_c (xi=0.6)')
plt.hlines([np.pi/4, -np.pi/4], -1.0, 1.0, lw=1, alpha=0.3)
plt.legend(loc='upper left');
plt.xlim(-0.002, 0.102)
plt.title('Bowshock shape parameters')
plt.xlabel('beta')
Note that all the shapes are hyperbolae for the CRW isotropic models, so $\theta_\mathrm{c} < 0$. For the proplyds, on the other hand, $\theta_\mathrm{c} > 0$ for small values of $\beta$.
Dmax = 15
mask = mm15
g = sns.jointplot(A(beta[mask]), np.degrees(thc(beta[mask], xi=0.8)),
stat_func=None, kind='hex', size=6,
xlim=(1.5, 3.0), ylim=(-90, 90))
g.set_axis_labels('A', 'th_c, degrees')
nexpect = mask.sum()/multiplier
g.ax_joint.text(1.55, 0.0,
"Sample with true separation $D < {}''$, $N = {:.1f}$".format(Dmax, nexpect),
fontsize='x-large');
Aprime, qprime = [], []
Aprime_p, qprime_p = [], []
for inc, b in zip(bigp.inc[mask], beta.value[mask]):
q = np.sqrt(b)/(1.0 + np.sqrt(b))
c = conproj_utils.Conic(A=A(b), th_conic=np.degrees(thc(b)))
Aprime.append(c.Aprime(inc))
qprime.append(q*c.g(inc))
cp = conproj_utils.Conic(A=A(b), th_conic=np.degrees(thc(b, xi=0.8)))
Aprime_p.append(cp.Aprime(inc))
qprime_p.append(q*cp.g(inc))
Aprime = np.array(Aprime)
qprime = np.array(qprime)
Aprime_p = np.array(Aprime_p)
qprime_p = np.array(qprime_p)
goodmask = np.isfinite(Aprime)
goodmask_p = np.isfinite(Aprime_p)
print(stats.describe(Aprime[goodmask])[:4])
print(stats.describe(Aprime_p[goodmask_p])[:4])
print(stats.describe(qprime[goodmask])[:4])
print(stats.describe(qprime_p[goodmask_p])[:4])
badmask = qprime > 1.0
badmask_p = qprime_p > 1.0
s = "Arc fraction = {:.3f}, q' > 1 fraction = {:.3f}"
print("Isotropic:", s.format(goodmask.sum()/len(goodmask),
badmask.sum()/len(goodmask)))
print("Proplyd:", s.format(goodmask_p.sum()/len(goodmask_p),
badmask_p.sum()/len(goodmask_p)))
sns.jointplot(beta.value[mask][goodmask_p],
np.abs(bigp.inc[mask][goodmask_p]),
kind='hex')
good = goodmask & ~badmask
good_p = goodmask_p & ~badmask_p
hist_kws = {"range": (1.0, 3.5), "alpha":0.6}
with sns.color_palette("coolwarm_r", 2):
ax = sns.distplot(
Aprime_p[good_p], label="proplyd",
kde=False, bins=50, hist_kws=hist_kws)
sns.distplot(
Aprime[good], label="isotropic",
ax=ax,
kde=False, bins=50, hist_kws=hist_kws)
ax.set_xlabel("A'")
ax.legend();
hist_kws = {"range": (0.0, 1.0), "alpha": 0.6}
with sns.color_palette("coolwarm_r", 2):
ax = sns.distplot(qprime_p[good_p],
label="proplyd",
kde=False, bins=50, hist_kws=hist_kws)
sns.distplot(qprime[good],
label="isotropic",
ax=ax, kde=False, bins=50, hist_kws=hist_kws)
ax.set_xlabel("q'")
ax.legend();
from functools import reduce
reduce?
Different ways of merging multiple dicts. They all work.
def dmerge(*dicts):
"""Merge an arbitrary number of dicts together.
Uses a simple `update` technique. Later dicts in
the argument list override earlier ones if there are
repeated keys.
"""
result = {}
for d in dicts:
result.update(d)
return result
def dmerge2(*dicts):
"See http://dietbuddha.blogspot.mx/2013/04/python-expression-idiom-merging.html"
return {k:v for d in dicts for k, v in d.items()}
def dmerge3(*dicts):
"See Power Baker's Feb 2014 comment on above page"
return reduce(lambda x, y: x.update(y) or x, dicts, {})
dmerge({"a": 1}, {"b": 2}, dict(c=3, d=4))
cmap_i = sns.blend_palette(["MidnightBlue", "PowderBlue"],
256, as_cmap=True)
cmap_p = sns.blend_palette(["IndianRed", "white"],
256, as_cmap=True)
nsample = 100
istart = 271
scatter_kws={
"alpha": 0.7,
"vmin": 0.0, "vmax": 0.9,
}
iso_kws={
"s": 2*bigp.rAU[mask][good][istart:istart+nsample],
"c": np.abs(np.sin(np.radians(bigp.inc[mask][good][istart:istart+nsample]))),
"cmap": cmap_i
}
prop_kws={
"s": 2*bigp.rAU[mask][good_p][istart:istart+nsample],
"c": np.abs(np.sin(np.radians(bigp.inc[mask][good_p][istart:istart+nsample]))),
"cmap": cmap_p
}
g = sns.jointplot(
qprime[good][istart:istart+nsample],
Aprime[good][istart:istart+nsample],
stat_func=None, kind='scatter', size=8,
xlim=(0.0, 0.5), ylim=(0.0, 3.5),
joint_kws=dmerge(scatter_kws, iso_kws))
g.ax_joint.scatter(
qprime_p[good_p][istart:istart+nsample],
Aprime_p[good_p][istart:istart+nsample],
**dmerge(scatter_kws, prop_kws)
)
sns.distplot(qprime_p[good_p][istart:istart+nsample],
kde=False, ax=g.ax_marg_x, color="IndianRed")
sns.distplot(Aprime_p[good_p][istart:istart+nsample],
kde=False, ax=g.ax_marg_y, vertical=True, color="IndianRed")
g.set_axis_labels('q\'', 'A\'')
rmin, rmax = stats.describe(bigp.rAU[mask][good_p][istart:istart+nsample])[1]
g.ax_joint.text(
0.01, 0.1,
"Sample with true separation $D < {}''$, $N = {:d}$, $r = {:.0f}$–${:.0f}$ AU".format(Dmax, nsample, rmin, rmax),
fontsize='x-large')
g.fig.savefig('monte-carlo-Aprime-qprime.pdf');
import matplotlib
matplotlib.colors.SymLogNorm?
plt.plot(conic.x(trange), conic.y(trange))
plt.axis([-3, 3, -3, 3])
Scratch area¶
sns.palplot(sns.color_palette("colorblind"))
sns.jointplot?
repr(jjdata[1]['Dprime_HA'])
a = dict(a=1, b=2)
b = {'x': 2}
c = a.copy()
c.update(b)
c
import astropy
astropy.__version__
from astropy import conf
conf.max_lines = -1
astropy.io.ascii.write?
(3e14*u.cm).to(u.au)
q = (1e14*u.cm)/(1.0*u.au)
"{0.value:0.03e} {0.unit:latex}".format(q)
plt.hist?
jjdata
data.pformat?
from astropy.io import ascii
ascii.write(data, Writer=ascii.Latex)
ascii.Latex?
ascii.latexdicts
data[data['ID'] == '177-341']['ID']
sns.jointplot?
sizedata[::10]
coord.SkyCoord.to_string?
coords0[:4].to_string('hmsdms')
coord.Angle.to_string?
coords0[:4].ra.to_string(unit=u.hourangle)
plt.scatter?
sns.jointplot?
import matplotlib
matplotlib.markers?
plt.colorbar?
f.axes[0] == fig._ax1
r = np.random.RandomState(5)
r.random_sample(10)
np.r_[[1, 2, 3], 0, 0, [1, 2]]
import seaborn
seaborn.palettes.dark_palette('red', as_cmap=True)
sns.palplot(sns.color_palette("RdPu", 5))
sns.color_palette("RdPu", 5)
sns.__version__
matplotlib.collections.Collection?